AD620003 


TECHNICAL  REPORT  ECOM-00392-1 


THE  INFLUENCE  OF  LOCAL  WINDS 

ON  FALLOUT 

SEMI-ANNUAL  PROGRESS  REPORT  NO.  1 


ECOif 


UNITED  STATES  ARMY  ELECTRONICS  COMMAND  •  FORT  MONMOUTH,  N.J 


CONTRACT  NO.  DA  28-043  AMC-00392(E) 

THE  TRAVELERS  RESEARCH  CENTER,  INC. 

250  CONSTITUTION  ^LAZA,  HARTFORD.  CONNECTICUT  06103 


Qualified  requesters  may  obtain  copies 
of  this  report  from  DDC . 

This  report  has  been  released  to  CFST1. 


Technical  Report  ECOM-00392-1 


August  1965 


THE  INFLUENCE  OF  LOCAL  WINDS 
ON  FALLOUT 


Semi-annual  Progress  Report  No.  1 
15  October  1964—15  April  1965 


Contract  No.  DA  28-043  AMC-00392(E) 
DA  Project  No.  1V0-14501-B-53A-O1 


Prepared  by 


James  D.  Kangos 
and 

Joseph  P.  Gerrity,  Jr. 


7464-164 


THE  TRAVELERS  RESEARCH  CENTER,  INC. 
250  Constitution  Plaza  Hartford,  Connecticut  06103 

for 


U.S.  ARMY  ELECTRONICS  COMMAND,  Fort  Monmouth,  N.  J. 


PURPOSE 


The  purpose  of  the  study  is  to  formulate  existing  information  on  the 
structure  of  local  wind  systems  and  local  (close-in)  fallout  characteristics 
into  a  computerized  computational  model  for  use  in  estimating  the  influence 
of  local  wind  circulations  on  the  deposition  of  close-in  fallout. 


ABSTRACT 


Information  on  the  structure  of  the  sea  breeze  obtained  by  previous  observa¬ 
tional  and  theoretical  studies  has  been  employed  to  formulate  a  kinematic  model 
of  the  wind  field  associated  with  the  sea  breeze.  Existing  information  on  the 
structure  of  the  stabilized  fallout  cloud  has  been  compiled  for  use  in  the  numerical 
computation  of  the  deposition  of  the  cloud  particles  through  specified  local  wind 
fields. 

A  program  was  constructed  for  the  numerical  computation  of  fallout  trans¬ 
port  through  an  analytically  prescribed  two-dimensional  local  wrind  circulation 
and  for  the  computation  of  its  ground  deposition  pattern.  The  computer  program 
was  used  to  evaluate  the  fallout  deposition  pattern  from  simple  models  of  the 
stabilized  cloud  through  alternative  models  of  the  sea-breeze  wind  fields. 

The  planned  development  of  a  three-dimensional  kinematic  model  of  the 
sea-breeze  wind  field  and  of  a  computer  program  for  the  computation  of  fallout 
transport  through  it,  is  discussed. 
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1.0  INTRODUCTION 


The  occurrence  of  a  nuclear  detonation,  either  in  the  air  or  under  the  sur¬ 
face,  with  subsequent  venting,  introduces  into  the  atmosphere  material  which  has 
been  rendered  radioactive-  This  material  is  distributed  through  the  atmosphere 
in  the  immediate  proximity  of  the  detonation  by  the  convective  circulation  induced 
through  the  heating  of  the  air  by  the  energy  released  in  the  detonation.  The 
extent  of  the  penetration  of  the  radioactive  material  through  the  atmosphere  is 
dependent  upon  the  energy  yield  of  the  detonation,  the  structure  of  the  atmosphere, 
the  nature  of  the  material  rendered  radioactive,  and  the  altitude  at  which  the 
detonation  occurred.  During  the  time  period  required  for  the  detonation- induced 
convection  to  subside,  complex  thermodynamic  and  radiochemical  processes  are 
active  within  the  mixture  of  air  and  radioactive  materials.  The  net  result  of 
this  period  of  radioactive  cloud  growth  is  the  distribution  of  the  particulate,  radio¬ 
active  material  throughout  the  region  occupied  by  the  cloud. 

In  the  past,  simplified  models  of  this  final  cloud  distribution  of  radio¬ 
active  material  have  been  employed  to  specify  the  initial  data  for  the  subsequent 
computation  of  the  redistribution  and  ultimate  deposition  of  the  particulate  matter 
upon  the  earth’s  surface.  At  the  present  time,  considerable  effort  is  being  expended 
to  improve  our  knowledge  of  the  mechanisms  relating  certain  of  the  factors  con¬ 
trolling  the  cloud  growth  to  the  characteristics  of  the  final  stage  in  the  process 
(widely  referred  to  as  the  “stabilized  fallout  cloud”). 

In  view  of  the  effort  to  improve  our  current  capability  to  specify  the  “source 
cloud”  characteristics,  it  is  considered  timely  to  improve  that  phase  of  the  fallout 
computation  which  transforms  the  “source  cloud”  data  into  subsequent  deposition 
patterns  or  transit  dose-rate  records.  One  aspect  of  this  transport  problem  which 
has  received  only  very  limited  consideration  thus  far  is  the  influence  of  the  so- 
called  “local  atmospheric  circulations”. 

The  work  being  reported  in  this  paper  is  an  extension  of  the  preliminary 
investigations  of  the  significance  of  local  circulations  (5,  6)  carried  out  at  The 
Travelers  Research  Center,  Inc.  (TRC),  under  contract  with  the  U.S.  Army 
Electronics  Laboratory. 
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This  earlier  work  led  to  certain  conclusions  which  we  will  now  briefly 
summarize; 

(a)  Two-dimensional,  vertical  circulations  which  do  not  change 
during  the  time  required  for  deposition  of  a  fallout  cloud  can  alter  the 
time  of  deposition  but  not  the  final  pattern  of  deposition. 

(b)  When  such  circulations  intensify  or  diminish  during  the  time 
required  for  deposition,  the  onset-time,  duration,  and  final  pattern  of  fall¬ 
out  deposition  can  be  modified. 

(c)  Two-dimensional,  horizontal  circulations  affect  the  horizontal 
pattern,  but  not  the  time  of  deposition.  Vertical  air  currents  must  exist 
if  the  time  of  deposition  is  to  be  modified. 

(d)  Three-dimensional  circulations,  even  when  temporally  steady, 
influence  both  the  time  and  pattern  of  deposition. 

(e)  The  compressibility  of  the  atmosphere  is  only  significant  when 
the  speed  of  the  vertical  air  currents  is  comparable  to  the  fall  speed 

of  the  fallout  particles. 

(f)  The  influence  of  the  two-dimensional,  temporally- varying, 
vertical  circulation  on  the  time  and  pattern  of  deposition  was  found  to 
decrease  by  an  order  of  magnitude  as  the  fallout  particle  fall  speed 
varied  from  25  cm  sec  1  to  1  m  sec  1 .  In  the  case  of  the  slower  falling 
particles,  significant  changes  (of  the  order  of  100%)  in  the  deposition 
pattern  and  time  of  deposition  (several  hours)  were  computed  over  limited 
regions. 

(g)  The  potential  significance  of  the  local  circulation  on  the 
pattern  and  time  of  deposition  depends  predominantly  on  the  time  and 
position  of  the  nuclear  detonation  in  relation  to  the  temporal  phase  in 
the  life-cycle  of  the  local  circulation. 

The  technical  objectives  of  the  current  program  are: 

(a)  <o  automate  the  procedures  for  computing  fallout  trajectories 
employee  in  the  past  elfort,  and  to  extend  the  method  so  that  general  three- 
dimensi  <nal  wind  circulations  may  be  accommodated;  and 


(b)  to  improve  upon  the  realism  of  the  two-dimensional  kinematic 
model  of  the  sea-breeze  circulation  utilized  in  the  previous  work;  to  extend 
the  scope  of  the  wind  models  treated  to  include  three-dimensional  sea- 
breeze  circulations  with  over-riding  synoptic -sc ale  flow  patterns,  and 
to  formulate  kinematic  models  of  the  mountain-valley  wind  field. 

In  this  progress  report,  we  will  present  the  automated  two-dimensional  sea- 
breeze  model  and  fallout  computation  programs.  The  results  of  a  limited  number 
of  numerical  computations  are  also  presented  and  discussed.  We  have  also  included 
discussions  of  the  three-dimensional  fallout  computation  logic  as  well  as  the  three- 
dimensional  sea-breeze  models  which  have  been  formulated  for  subsequent  use. 

In  Appendix  B,  we  present  the  computer  program  which  has  been  developed  for 
computing  the  deposition  of  fallout  passing  through  a  two-dimensional  local 
circulation. 
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2,0  THE  SEA  BREEZE 


2.1  Review  of  Literature 

Up  to  the  late  1950’s,  observations  of  the  sea  breeze  were  confined  to 
vertical  soundings  of  the  horizontal  wind  components  at  a  single  geographic 
location,  and  to  routine  surface  observations  of  meteorological  elements.  A  brief 
summary  of  these  observations  based  on  an  article  by  Defant  (2]  follows. 

When  the  prevailing  large-scale  flow  is  light,  the  sea  breeze  develops  as  a. 
small  circulation  in  the  immediate  vicinity  of  the  coast.  The  circulation  then 
gradually  increases  in  depth  and  horizontal  extent,  both  landward  and  seaward. 

The  existence  of  the  large-scale  synoptic  flow  from  land  to  sea  results  in  an  initial 
development  out  at  sea.  The  intensification  and  advance  of  the  sea-  breeze  land¬ 
ward  is  much  more  gradual  and  reaches  the  coast  later  in  the  afternoon,  similar 
to  a  front  with  a  characteristic  wind  shift. 

The  sea  breeze  consists  of  a  landward  current  adjacent  to  the  earth’s  surface 
and  a  much  weaker,  but  deeper,  return  flow  above.  The  top  of  the  landward  current 
ranges  from  a  few  hundred  meters  on  a  mid-latitude  coast  to  approximately  2  km 
along  tropical  coasts.  The  landward  penetration  of  the  sea  breeze  depends  on  the 
land— sea  temperature  difference  and  the  large-scale  synoptic  flow.  During  ideal 
conditions  (weak  large-scale  flow),,  the  sea  breeze  is  frequently  observed  as  far 
inland  as  30  to  50  km  in  mid-latitude  regions  and  sometimes  several  hundred 
kilometers  inland  in  tropical  regions. 

The  maximum  horizontal  w'ind  velocities  in  a  sea-breeze  circulation  are 
-1 

approximately  10  m  sec  , whereas  the  vertical  velocities  are  about  tw'O  orders  of 
magnitude  smaller.  The  daily  maximum  wmds  generally  occur  shortly  after  the 
temperature  maximum. 

During  the  latter  portion  of  the  1950’s,  Fisher  [7,9]  conducted  a  series  of 
observational  studies  of  the  sea  breeze  on  the  coast  of  New  York  and  Rhode  Island. 
For  the  first  time,  a  simultaneous  set  of  vertical  soundings  of  the  horizontal  wind 
and  temperature  were  obtained  during  the  presence  of  a  sea  breeze.  A  series  of 
observation  sites,  aligned  along  a  line  normal  to  the  coast,  provided  the  data  net¬ 
work  from  w'hich  a  vertical  cross-section  of  the  structure  of  the  sea  breeze  could 
be  constructed. 
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Fisher’s  observational  program,  which  was  sponsored  by  the  U.S,  Army 
Signal  Research  and  Development  Laboratory,  Fort  Monmouth,  New  Jersey,  in 
addition  to  substantiating  many  previously  known  characteristics  of  the  sea- 
breeze  (  'tculation,  revealed  certain  characteristics  which  had  not  been  noted 
before.  These  new  facets  of  the  sea  breeze  had  escaped  previous  theoretical 
investigations  based  upon  linear  analytical  mathematics. 

Both  Fisher  [8,9]  and  Estoque  [3,4]  formulated  physical  models  of  the  sea 
breeze  for  study  by  means  of  numerical  integration  using  electronic  computers. 
Fisher’s  numerical  experiments  succeeded  in  reproducing  the  principal  character¬ 
istics  of  a  sea  breeze  in  a  calm  environment.  They  also  reproduced  the  tendency 
for  a  sudden  acceleration  of  the  onshore  flow  that  was  observed  to  occur  during 
the  declining  stage  of  the  circulation’s  life  cycle. 

Estoque’s  model  was  developed  somewhat  differently  than  Fisher’s,  with 
the  intention  of  investigating  the  influence  of  variations  in  the  over-riding  synoptic 
flow  pattern  on  the  development  of  the  sea-breeze  circulation.  The  most  sig¬ 
nificant  characteristic  of  the  circulation  found  by  Estoque  is  the  sea-breeze  front. 
This  sharp,  leading  edge  of  the  sea  breeze  was  computed  when  the  overriding 
synoptic  flow  was  set-up  to  oppose  the  landward  penetration  of  the  sea  breeze. 

The  vertical  component  of  the  velocity  field  was  a  maximum  in  these  cases. 

Its  magnitude  was  as  large  as  20  cm  sec  *. 

Although  the  numerical  results  obtained  by  Estoque  and  Fisher  are  most 
interesting,  they  do  not  provide  a  simple  representation  of  the  sea-breeze  circula¬ 
tion  in  terms  of  observable  parameters. 

2.2  Two-dimensional  Model 

In  the  first  half  of  this  study,  a  two-dimensional  sea-breeze  model  was  used 
in  an  attempt  to  determine  if  the  results  of  the  proposed  fallout  computation 
scheme  were  reasonable.  This  model  was  developed  by  utilizing  all  available 
observational  material  [7,  9]  in  an  effort  to  simulate  the  sea  breeze  as  closely  as 
possible  in  an  x,  z-plane.  To  take  into  account  the  development  of  the  sea  breeze, 
the  parameters  influencing  the  wind  field  are  considered  functions  of  time. 

Because  parameters  governing  the  position  and  intensity  of  the  sea  brcc  zc  3.rc  time- 
dependent  in  our  two-dimensional  model,  the  observational  data  of  Fisher  [7,9] 


were  carefully  studied  so  that  their  behavior  could  be  expressed  by  simple 
functions.  The  parameters  important  in  the  development  are: 

(a)  X(t),  the  location  of  the  vertical  center  of  the  circulation  in 
relation  to  the  coast, 

(b)  ZT(t),  height  of  the  top  of  the  circulation, 

(c)  U(t),  maximum  value  of  the  u -component  of  the  wind  in  the  local 
circulation,  and 

(d)  cr( t),  the  horizontal  extent  of  the  sea  breeze. 

Three  parameters  which  appear  in  our  basic  equation  as  functions  of  time  will  be 
treated  as  constants  as  they  are  nearly  constant  according  to  observations.  These 
are: 

(a)  Z^(t),  the  height  of  the  base  of  the  local  circulation, 

(b)  ^^(t),  height  of  the  maximum  value  of  the  u-component  in  the 
local  circulation,  and 

(c)  Zp(t),  the  height  of  the  top  of  the  layer  influenced  by  the  sea 
breeze. 

The  values  of  these  parameters  are  0.07  km,  0.4  km,  and  3.0  km  respectively. 
A  set  of  linear  and  quadratic  functions  was  developed  approximating  the  behavior 
of  the  other  parameters.  Employing  these  functions,  the  development  of  the  sea 
breeze  begins  at  hour  1000  local  time,  and  intensifies  and  remains  on  the  coast 
until  hour  1200.  The  sea  breeze  then  begins  to  progress  inland  and  continues 
intensification  until  hour  1500  when  it  has  reached  peak  development  at  a  distance 
of  12  km  inland.  From  1500  to  1800  hours,  the  sea  breeze  dissipates  and  retro¬ 
gresses  toward  the  coast.  At  1800,  it  is  again  centered  on  the  coast,  where  it 
remains  until  it  is  completely  dissipated  at  2000.  Incorporating  these  functions 
into  the  model  provides  us  with  a  good  approximation  of  a  sea  breeze  during  its 
entire  life  cycle.  As  it  is  impossible  to  specify  the  behavior  of  X(t)  and  Z^,(t)  with 
a  single  equation,  the  expressions  are  presented  with  the  time  periods  in  which  they 
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are  applicable.  For  convenience  in  evaluating  the  results,  the  coordinate  system 
is  centered  on  the  coast.  These  functions  are  as  follows; 


(a)  X,  position  of  center  of  circulation. 
10  £  t  <  12 
18  <  t  20 


X(t)  =  0 


X(t)  =  (~  1.33t2  +  40t  -  288)  105, 


12  *  t  as  18. 


(b)  Z^t),  height  of  the  top  of  the  local  circulation. 


ZT(t)  =  (0.4t  -  4)10  , 


10  s  t  <  12 


ZT(t)  ®  (-  +  1.33t  -  8.8)  105, 


12 


18 


ZT(t)  =  {-0.4t  +  8)10  ,  18  <  t  ^  20. 

(c)  U(t),  u-component  of  the  wind  in  the  local  circulation. 


U(t)  =  (-  ~  +  4.1  t  -  22.6)  102,  10  S  ts  20. 

(d)  cr(t),  horizontal  extent  of  the  sea  breeze. 

ff(t)  =  (-0.88t2  +  26. 4t  -  166)  IQ5,  10  £  t  ^  20. 


(2-1) 


(2-2) 


(2-3) 

(2-4) 

(2-5) 

(2-6) 

(2-7) 


By  using  these  expressions  in  the  sea-breeze  model  presented  below,  wind 
fields  approximating  a  two-dimensional  sea  breeze  can  be  obtained.  Figures  2-1 
and  2-2  illustrate  the  sea  breeze  circulation  for  hours  1100  and  1500  respectively. 

Because  it  is  often  found  that  an  over-riding  synoptic  flow  is  present  during 
the  development  of  a  sea  breeze,  the  model  presented  will  contain  an  undisturbed 
component  of  the  wind,  and  a  weak  return  flow  above  the  local  circulation. 

Now,  introducing  a  simple  flow  in  a  z,t-plane,  we  consider  a  wind  field 
composed  of  an  undisturbed  part  and  a  local  circulation.  The  undisturbed  part 
is  assumed  to  be  solely  a  component  of  velocity  in  the  x-direction.  The  local 
circulation  is  assumed  to  obey  the  equation. 


dw  da 

dz  9x ' 


(2-8) 
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Height, 


Distance,  km 

Fig.  2-2.  u  and  w  components  of  sea  breeze  at  time  1500  {u  (m  sec  *)>  w  (cm  sec  )]. 
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In  this  case,  we  will  specify  the  u-velocity  fi^ld  and  use  it  and  an  arbitrary 
boundary  condition,  to  solve  w  from  Eq.  (2-8), 


w  -  0  at  z  =  z.  for  all  x,t. 
h 


(2-9) 


At  this  point,  it  is  convenient  to  introduce  the  notation  4»(x,t),  which  is  a  weight¬ 
ing  function  whose  value  ranges  between  0  and  1: 

.2  1 


4»(x,t)  *  exp 


{- 


d&(x,  t)  _ 


{ 


2{o(t)}2  / 

x  -  X(t)l  *(x,t). 

«t>2  > 


(2-10) 

(2-11) 


The  total  velocity  components  will  be  denoted  as  u  and  w.  The  undisturbed 
x-component  will  be  denoted  by  u,  the  local  wind  x -component  by  u,  and  the  local 
wind  z -component  by  w. 

The  equation  for  the  undisturbed  flow  is; 

e-<az>i 


u(z)  *  G 


<-[ 


J 


cos  (oz)}. 


where  G  is  the  geostrcphic  wind  velocity,  and  a  is  the  Ekman  parameter 
/  1/2 

(f/ 2k)  .  The  equations  governing  the  local  circulation  are: 


fi(x,z,t)  »  0, 


Z  ^  Zh(t) 


fi(x.z.t)  * 


U(t)*<x,t)  (Z  -  Zh(t)J 

~  izm(t)  -  zh(t)] 


Z  <  Z  ^  Z  (t) 
n  m 


(2-12) 


(2-13) 

(2-14) 


u(x,z,t)  - 


U(t)*(x,t)  tZT(t)  -  Z] 

lZT(t)  -  zm<t)j2 


d(x,z,t)  *  -  N(x,t)  F(x,t), 


Zm(t)  <  Z  £  ZT(t)  (2-15) 


ZT(t)  <  Z  s  Zp(t)  (2-16) 


where  N(x,t)  represents  the  net  flow  in  the  positive  x-direction  between  Z,  and 

h 

Z^,,  and  is  given  by 


N(x,t)  =  U(t)  *(x 


A 


Zm(t>  zh<‘>  M'-yi], 


(2-17) 


10 


The  function,  F(z,t),  distributes  the  compensating  return  flow  normally  through* 

out  the  depth  Z_  to  Z  ,  and  is  given  by 
I  r 


in  which 


F(z*t)  *  exp 


o(t)  » 


Zp(t)  -  zT<t) 


IZ  -  Z(t)] 
2 

2cr(t) 


(2-18) 


and 


Z(t)  - 


ZF(t)  +  ZT(t) 


Now  from  Eq.  (2-8)  and  Eqs.  (2-12)  through  (2-16)  we  find  that 


and 


9w 

9u  = 

dz 

dx 

dw 

du  _ 

dz 

9x 

Z  s  Zh(t) 


o(t)z 


u(x,z,t),  Zh(t)  <  Z  s  Zp(t) 


(2-19) 

(2-20) 


Using  the  boundary  condition  w  =  0  at  Z  -  Z^(t),  we  now  integrate  Eq.  (2-20)  to 

obtain  w  in  the  region  between  Z,  and  Z  . 

n  r 


Writing 


w  =  w  [Z  =  Z  (t)j  , 
m  mv  ' 


we  find 


w(x,z)  = 


X  -  X(t) 
2 

a(t) 


j  *(x. 


t)  U(t) 


[Z  -  Z  (t)J 
_ m 

2[Zm(t)  -  Zh(t)]* 


w(x,  Z)  =  ft 


m 


Zh(t)  <  Z  s  Zffi(t) 

X  -  Jjt)  _ *tX»t).U(t)  _ 

o(t)2  3[ZT(t)  -  Zm(t)]2 


[ZT(t)  -  Z]' 


lZT(t)  ' 


Z  (t)  <  Z  ^  Z  (t). 
m  l 


(2-21) 


(2-22) 


Now  writing 


ft  =  ft  [Z  =  Z_(t)J, 


we  find  . 

,  -  N(x.z)  i 

W(x,  Z)  =  WT  -  exp 


Z 


[- 


jZ  -  Z(t)l 
21a(t)]2 


2  i 


dz 


X  -  %(t) 


]• 


<7(t) 

ZT(t)  <  Z  s  ZF(t). 


(2-23) 
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t  Three-dimensional  Model 

As  was  mentioned  earlier,  local  circulations  are  three-dimensional  in 
nature  and  their  dimensions  are  strongly  time  dependent.  The  two-dimensional 
model  described  above  was  developed  to  test  the  computation  scheme  used.  How¬ 
ever,  our  final  objective  in  this  study  is  to  develop  three-dimensional,  kinematic 
local  circulation  models  and  to  compute  the  influence  of  such  local  winds  on  the  pattern 
of  close-in  fallout  deposition.  This  requires  that  we  be  able  to  describe  quantita¬ 
tively  the  three  velocity  components  of  natural  antitriptic  wind  fields  in  terms  of 
a  limited  number  of  observable  parameters.  Because  observational  material  of 
the  sea  breeze  in  three  dimensions  are  non-existent,  we  have  turned  to  the  results 
of  theoretical  studies  in  an  effort  to  seek  such  representations  of  the  sea  breeze. 

In  the  work  of  Pierson  [14]  and  Haurwitz  [11],  we  have  found  compatible  studies 
of  temporal  and  spatial  structure  of  the  sea  breeze.  Although  the  formulas  derived 
from  their  work  are  complicated,  we  must  emphasize  that  they  do  not  admit 
representation  of  the  non-linear  characteristics  of  the  sea  breeze  found  in  the 
physical  models  of  Estoque  and  Fisher.  The  generation  of  three-dimensional,  time- 
dependent  wind  fields  from  the  theoretical  solutions  obtained  both  analytically  and 
numerically  will  require  the  imaginative  manipulation  of  the  solutions  to  fit  the 
limited  observational  data.  Although  we  have  initiated  this  work,  it  would  be  pre¬ 
mature  to  report  upon  it  in  detail  at  this  time. 
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3.0  THE  FALLOUT  CLOUD 


3.1  Review  of  Literature 

Since  the  advent  of  atomic  weapons,  numerous  studies  have  been  conducted 
dealing  with  the  problem  of  both  close-in  and  distant  fallout.  To  date,  one  of 
the  most  perplexing  problems  still  confronting  investigators  is  the  distribution 
of  radioactive  particles  within  the  stem  and  mushroom  cloud.  Kellogg,  Rapp,  and 
Greenfield  112]  made  some  estimates  of  particle  distribution  based  on  observed 
fallout  and  a  reconstruction  of  what  the  initial  distribution  must  have  been.  How¬ 
ever,  in  reconstructing  the  initial  distribution,  two  assumptions  were  made: 

(a)  Ten  percent  of  the  debris  was  distributed  evenly  throughout 
the  stem  and  the  remaining  portion  was  in  the  mushroom  cloud. 

(b)  The  particles  in  the  cloud  were  distributed  in  such  a  way  that 
there  was  a  constant  mixing  ratio  between  the  particles  and  the  air. 

This  would  mean  that  there  was  complete  mixing  in  the  mushroom  cloud  of  90 
percent  of  the  particles  and,  therefore,  the  density  of  the  particles  drops  off 
exponentially  with  height  in  exactly  the  same  way  as  the  air  density,  in  the  stem 
this  is  not  true,  and  it  can  be  seen  that  a  constant  density  of  particles  with  height 
implies  a  mixing  ratio  that  increases  with  height.  That  is,  relative  to  the  air, 
there  are  more  particles  in  the  top  of  the  stem  than  in  the  lower  part. 

Anderson  (1  ],  in  a  theoretical  study,  found  that  close-in  fallout  is  comprised 
mainly  of  particles  with  significant  gravitational  fall  velocities  which  are  best 
defined  in  terms  of  a  minimum  particle  size  of  about  50  microns  in  diameter. 
Particles  smaller  than  this  exhibit  a  rapid  decrease  of  fall  speed  with  decrease 
of  size.  Consequently,  the  particle  sizes  important  in  close-in  fallout  range 
from  500  microns  to  25  microns  in  radius. 

Machta  {13],  in  a  study  of  the  dimensions  of  an  atomic  cloud,  developed  a  semi- 
empirical  theory  relating  the  dimensions  to  atmospheric  stability  and  the  rate  of 
entrainment.  This  semi-empirical  theory,  although  widely  used,  provides  only 
rough  estimates;  however,  as  there  is  a  current  lack  of  a  more  accurate  theory, 
cloud  dimensions  used  in  our  study  are  based  on  it. 
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3.2  Computation  of  Fallout  in  a  Two-dimensional  Model 

In  computing  the  motion  of  a  fallout  cloud  through  a  sea-breeze  circulation, 
the  governing  factors  are  the  horizontal  and  vertical  wind  velocities  and  the  fall 
speed  of  the  particles  contained  within  the  cloud.  In  formulating  this  portion  of 
the  problem,  we  considered  a  rather  simple  scheme  whereby  the  fallout  cloud 
was  divided  into  n  slices.  Each  slice  contained  a  specified  fraction  of  the  total 
mass  which  was  assumed  to  be  uniformally  distributed  throughout  the  slice.  At 
this  stage,  we  prescribed  a  finite  number  of  points  describing  the  boundary  of 
each  slice.  The  slices  are  then  allowed  to  descend  through  the  local  circulation 
with  a  given  fall  velocity  until  all  the  slices  have  been  deposited.  This  procedure 
is  repeated  using  a  spectrum  of  fall  velocities. 

In  determining  the  trajectory  of  the  cloud,  or  slices  in  our  case,  the  u  and  w 
wind  components  must  be  computed  at  the  points  describing  the  slices.  Having 
obtained  the  wind  components  at  these  points,  the  horizontal  and  vertical  displace¬ 
ments  of  the  points  are  computed  for  a  given  time  interval  (At)  by  using  the  follow¬ 
ing  formulas: 

Ax  =  uAt 

and 

Az  =  (w  +  VF)  At 

where  VF  =  fall  velocity  of  particles  (VF  <  0). 

After  computing  the  displacements  and  obtaining  a  new  set  of  coordinates 
for  the  points  describing  the  slice,  the  process  is  repeated  until  all  slices  have 
been  deposited  on  the  surface.  However,  it  is  possible  that  the  distance  between 
two  adjacent  points  on  the  slice  boundary  may  become  too  large,  thus  making  the 
results  unreliable.  Therefore,  a  set  of  limits  was  established  whereby  if  the 
distance  between  any  two  adjacent  points  exceeds  these  values,  the  iteration,  for 
that  particular  time  step  will  be  repeated  using  half  the  original  value  of  At.  If 
the  points  are  still  too  widely  separated  after  three  iterations,  the  program  is 
terminated  for  that  particular  fall  velocity. 

During  the  time  of  descent  through  the  local  circulation,  the  area  of  the  slice 
must  be  computed  after  each  time  step  because  the  area  may  change  due  to  con¬ 
traction  or  expansion,  depending  upon  the  displacement  of  the  points.  Having 
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specified  the  mass  within  the  slice,  the  density  can  be  computed.  At  the  time  deposi¬ 
tion  commences,  the  area  of  the  slice  found  below  the  ground  must  be  computed  so 
that  the  mass  of  material  deposited  can  be  determined  using  the  density  value  from 
the  previous  time  step.  The  fallout  is  then  assigned  to  that  particular  interval  on 
the  x-axis.  The  total  mass  is  now  modified  by  subtracting  the  amount  of  mass  that 
has  been  deposited.  The  area  ot  the  slice  remaining  above  the  ground  is  computed 
and,  using  the  modified  mass,  a  new  density  is  determined.  At  this  stage,  the  points 
that  descended  through  the  surface  are  now  replaced  by  points  on  the  surface 
vertically  above  them,  so  that  the  number  of  points  defining  the  slice  remains 
unchanged.  The  remaining  portion  of  the  initial  slice  is  now  resting  on  the  ground. 
The  u  and  w  wind  components  are  again  computed  for  all  the  points.  Once 
deposition  begins,  the  points  resting  on  the  surface  can  only  move  vertically  down¬ 
ward,  This  procedure  is  continued  until  the  slice  has  been  completely  deposited 
(see  Appendix  A). 

Because  the  fallout  cloud  was  initially  divided  into  n  slices,  each  slice 
must  be  treated  individually  until  all  slices  have  been  deposited. 

3.3  Computation  of  Fallout  in  the  Three-dimensional  Model 

In  formulating  the  three-dimensional  problem,  we  initially  considered  a 
rather  complex  scheme  for  computing  both  the  motion  of  the  fallout  cloud  and 
the  changes  in  the  concentration  of  the  fallout  within  the  cloud.  The  procedure 
would  generalize  the  two-dimensional  scheme  utilized  in  the  computations  reported 
above.  The  surface  of  the  fallout  cloud  would  be  approximated  by  a  polyhedron 
with  triangular  faces  and  the  vertices  of  the  polyhedron  cloud  would  be  associated 
with  several  of  the  triangular  faces.  By  proper  indexing  of  the  vertices  and  facets, 
the  interior  of  the  cloud  could  be  determined  throughout  the  trajectory  computa¬ 
tion.  A  principal  advantage  of  this  scheme  would  have  been  its  capacity  for 
admitting  computation  of  changes  in  the  fallout  concentration  within  the  cloud. 

Subsequent  consideration,  however,  led  us  to  decide  on  utilizing  the  simplifying 
assumption  that  the  fallout  concentration  was  conserved  during  the  motion  of  the 
cloud.  Consistent  with  thi3  assumption  is  a  simpler  description  of  the  cloud  and 
a  simpler  trajectory  computation  procedure,  somewhat  similar  to  that  used  in 
the  Ford  model  [15], 
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The  stabilized  cloud,  which  we  will  assume  to  be  eyclindricai  with  stem 
and  mushroom  cap,  is  divided  into  discs.  Each  disc  is  assumed  to  enclose  a 
specified  fraction  of  the  total  radioactive  mass  of  the  entire  cloud.  This  frac¬ 
tional  mass  content  is  further  assumed  to  be  distributed  according  to  a  particle 
size  spectrum  appropriate  to  the  position  of  the  disc  within  the  cloud. 

At  this  point,  a  further  assumption  is  made  that  the  continuous  particle 
size  spectrum  may  be  replaced  by  a  discrete  and  finite  line  spectrum.  A  unique 
fall  velocity  corresponds  to  each  line  in  the  spectrum.  Thus,  to  compute  the 
motion  of  the  material  originally  within  a  particular  disc,  the  disc  must  be  replaced 
by  as  many  sub-discs  (wafers)  as  there  are  lines  in  the  discrete  spectrum. 

When  the  center  of  gravity  of  the  disc  is  at  an  altitud'  above  which  local 
circulations  are  allowed,  we  move  each  wafer  with  the  wind  field  existing  at 
its  center  of  gravity.  When  the  center  of  a  wafer  reaches  an  altitude  below 
which  local  circulations  are  permitted,  we  further  subdivide*  the  wafer  into 
smaller  volume  elements.  Then  each  element  is  followed  through  the  local  circu¬ 
lation  until  deposition  occurs.  The  position  of  deposition  and  the  quantity  of 
material  deposited  is  recorded  and  the  computation  continued. 

*Refer  to  Appendix  C  for  subdivision  scheme. 
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4,0  COMPUTED  EFFECTS  OF  THE  SEA-BREEZE  CIRCULATION  ON  THE  FALL- 


OUT  DISTRIBUTION 
4.1  Background 

The  fallout  distribution  depends  upon  the  circulation  through  which  the 
particles  must  descend.  Several  experiments  were  performed  whereby  clouds  of 
various  dimensions  and  containing  the  spectrum  of  particle  sizes  were  simulated 
to  descend  through  the  circulation  produced  by  the  sea-breeze  model.  It  will  be 
noted  at  this  point  that,  in  these  experiments,  we  are  dealing  solely  with  the  mush¬ 
room  cloud  and  not  the  stem,  as  information  pertaining  to  its  dimensions  and 
particle  size  distribution  is  lacking.  The  dimensions  oi  the  mushroom  cloud, 
particle  size  distributions,  and  associated  terminal  velocities  were  obtained  from 
the  investigation  of  Kellogg,  Rapp,  and  Greenfield  (12).  Two  different  cloud  sizes 
were  used  in  the  experiments,  resulting  from  yields  of  10  kt  and  300  kt.  In  the 
10  kt  yield,  the  cloud  is  2.4  km  in  diameter  and  1.2  km  thick.  For  the  300  kt  yield, 
the  cloud  is  15  km  in  diameter  and  5  km  thick.  The  clouds  were  divided  into  discs 
0.61  km  and  1  km  thick,  respectively.  Therefore,  in  the  10  kt  and  300  kt  cases,  we 
computed  the  descent  of  two  and  five  discs  r-'  spectively,  through  the  local  circula¬ 
tion.  The  particle  sizes  used  in  both  cases  ranged  from  300  to  25  microns,  with 

-1  *1 

associated  terminal  velocities  of  500  cm  sec  to  25  cm  sec  *  The  base  of  the  clouds 
were  initially  at  3  km  above  the  surface . 

The  initial  position  of  the  cloud  in  relation  to  local  circulation  will  produce 
different  fallout  distributions,  and  the  experiments  were  designed  with  this  in 
mind.  Basically,  two  different  circulations  were  used,  (1)  sea-breeze  circulation 
imbedded  in  an  over-riding  synoptic  flow  from  the  sea  to  land,  and  (2)  a  sea- 
breeze  circulation  in  the  absence  of  an  over-riding  synoptic  flow.  The  positions  of 
the  10  kt  yield  experiments  were:  (1)  1.6  km  out  at  sea,  and  (2)  directly  over  the 
coast  line.  Similarly,  the  positions  of  the  300  kt  yield  experiments  were:  (1)  15  km 
out  at  sea,  and  (2)  directly  over  the  coast  line.  In  the  analysis  of  the  results,  the 
fallout  distributions  were  identical  in  both  cases  (cloud  located  offshore  and  on- 
coast),  except  that  they  were  shifted.  The  position  and  configuration  of  the  cloud 
is  a  variable  in  this  model,  (x-  and  z-  coordinates  describing  cloud)  and,  there¬ 
fore,  any  number  of  experiments  may  be  performed  with  the  computer  program. 
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To  investigate  the  effect  of  the  circulation  on  the  fallout  distribution,  the 
assumption  was  made  that  each  disc  was  a  homogeneous  layer  of  particles  with 
uniform  fall  velocities.  Nine  different  fall  velocities  were  used  and  are  presented 
in  Table  4-1. 


TABLE  4-1 
FALL  VELOCITIES 


Radius  (m) 

T'  rminal  velocity 
(cm  sec-1) 

300 

500 

250 

400 

175 

300 

125 

200 

90 

150 

70 

100 

60 

75 

40 

50 

25 

25 

Each  disc  was  allowed  to  descend  through  the  local  circulation  nine  times. 
Because  information  relating  the  mass  of  particles  contained  within  a  cloud  to 
weapon  yield  is  lacking,  a  value  of  100  was  used.  This  permits  us  to  express  the 
distribution  of  fallout  on  a  surface  as  a  percentage  of  the  total  mass  contained 
within  the  cloud. 

4.2  Local  Circulation  in  the  Absence  of  an  Overriding  Flow 

Two  experiments  were  performed  in  which  clouds  from  a  10  kt  and  300  kt 
yield  were  simulated  to  descend  through  a  sea-breeze  circulation  in  the  absence  of 
an  overriding  flow.  Figure  4-1  illustrates  the  distribution  of  fallout  as  a  result 
of  a  cloud  produced  from  a  10  kt  yield,  located  1.6  km  out  to  sea,  and  simulated  to 
descend  through  the  local  circulation  (solid  line)  and  the  resulting  distribution  in 
the  absence  of  any  circulation  (stipple  area).  In  this  case,  the  simulated  local 
circulation  modified  the  fallout  pattern  to  the  extent  that  approximately  50%  of  the 


18 


Absence  of  wind 


19 


Fig.  4-1.  Fallout  distribution  from  a  10-kt-yield  burst  located  offshore  with  no  overriding  flow. 


mass  was  distributed  outside  of  the  area  directly  beneath  the  initial  position  of 
the  cloud.  In  the  second  case  (see  Fig.  4-2),  a  cloud  from  a  300  kt  yield  and 
located  15  km  out  to  sea,  had  only  10%  of  its  mass  transported  away  from  the 
initial  position.  Particles  with  fall  velocities  of  less  than  100  cm  sec  1  were 
influenced  to  a  greater  degree  than  larger  particles,  as  was  shown  by  Feteris  [6], 
The  experiments  conducted  included  particles  with  fall  velocities  of  25  cm  sec  1 
However,  after  3  hours  of  real  time,  deposition  had  not  occurred  and  the  points 
described  in  the  boundary  of  the  slice  become  too  widely  separated  to  fall  within 
the  specified  limits  after  three  iterations;  therefore,  the  prograr.  was  terminated 
for  that  particular  fall  velocity.  This  happened  because  the  particles  with  fall 
velocities  of  25  cm  sec  1  and  less  are  greatly  influenced  by  the  local  circulation 
and,  had  we  assigned  extremely  large  limits,  deposition  would  have  occurred. 

The  results  however ,  would  have  been  unreliable . 

4,3  Local  Circulation  in  an  Overriding  Flow 

Two  experiments  were  also  conducted  using  a  local  circulation  imbedded  in 
an  overriding  flow.  The  positions  of  the  clouds  were,  as  before,  (1)  10  kt  yield, 
1.6  km  offshore,  and  (2)  300  kt  yield,  15  km  offshore. 

In  the  10  kt  yield  experiment,  the  fallout  distribution  was  quite  different 
from  the  pattern  obtained  from  the  local  circulation  only.  The  fallout  was  dis¬ 
tributed  over  a  larger  region  (27  km)  primarily  due  to  the  overriding  flow.  The 
distribution  of  fallout  was  also  computed  in  this  case  for  an  overriding  flow  with¬ 
out  the  local  circulation.  Figure  4-3  shows  that  the  local  circulation  did  influence 
the  fallout  pattern,  producing  peaks  and  minimums  in  relation  to  the  distribution 
obtained  solely  from  the  overriding  flow.  The  sea  breeze  in  this  case  did 
influence  the  fallout  pattern;  however,  it  should  be  kept  in  mind  that  with  a  very 
weak  overriding  flow  (of  the  order  of  1  cm  sec  \  the  local  circulation  may  have 
a  greater  influence  on  the  fallout  distribution. 

In  the  case  of  the  300  kt  yield  burst  (cloud  location  was  15  km  offshore), 
the  fallout  distribution  was  greatly  influenced  by  the  overriding  flow  (5  m  sec  1) 
because  the  area  covered  by  fallout  was  quite  large  (92  km).  Figure  4-4  shows 
that  a  marked  peak  exists  between  the  coast  and  the  initial  position  of  the  cloud. 
This  results  from  particle  sizes  with  fall  velocities  greater  than  100  cm  sec  \ 
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Fig.  4-2.  Fallout  distribution  from  a  300-kt-yield  burst  located  offshore  with  no  overriding  flow. 
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Fig.  4-3.  Fallout  distribution  from  a  10-kt-yield  burst  located  offshore  with  an  overriding  flow  of  5  m  sec 
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ig.  4-4.  Fallout  distribution  from  a  300-kt-yield  burst  located  offshore  with  an  overriding  flow  of  5  m  sec 


5.0  CONCLUSIONS 


The  primary  objective  in  this  study  is  to  develop  three-dimensional  local 
circulation  models  for  computation  of  close-in  fallout  patterns.  In  the  develop¬ 
ment  of  a  two-dimensional  sea-breeze  model,  it  has  been  demonstrated  that  such 
models  can  be  expressed  in  terms  of  analytical  functions  and  can  be  readily  pro¬ 
grammed  for  rapid  computation  of  fallout  distributions  as  influenced  by  such  cir¬ 
culations  .  Also,  the  results  obtained  by  Feteris  [6]  by  hand  computations  have  been 
verified  by  using  an  automated  two-dimensional  sea-breeze  model. 
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7.0  PROGRAM  FOR  THE  NEXT  INTERVAL 


As  indicated  above,  we  have  initiated  the  formulation  of  a  three-dimensional, 
time-dependent  kinematic  model  of  the  sea-breeze  wind  field.  This  work  will 
continue  and  very  likely  be  extended  to  consider  models  of  the  mountain- 
valley  wind.  A  new  computer  program,  based  on  the  logic  outlined  in  Section  3.3 
and  Appendix  C,  will  be  written.  Finally,  numerical  results  will  be  obtained  from 
experiments  judiciously  designed  to  indicate  the  influence  on  the  deposition  of  close- 
in  fallout  exerted  by  the  kinematic,  local  wind  models. 
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8.0  IDENTIFICATION  OF  PERSONNEL 


8,1  Extent  of  Participation 


TABLE  8-1 

PERSONNEL  PARTICIPATION* 


Name 

Title 

Total  hours  worked 
during  the  quarter 
(approximate) 

«J.  P.  Gerrity 

Principal  Investigator 

232 

J.  D.  Kangos 

Project  Scientist 

673 

E.  A.  Newburg 

Mathematician 

45 

P.  S.  Brown 

Mathematician 

38 

A.  F.  Saunders 

Programmer 

12 

J.  A.  Sekorski 

Programmer 

237 

♦Secretarial,  administrative,  and  drafting  assistance  was  also  provided. 


8.2  Biographies  of  Key  Personnel 
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Cwb-10759,  The  Travelers  Research  Center,  Hartford,  Nov.  1964. 


32 


Alexander  F.  Saunders,  B.S, 


Professional  Experience 

Mr.  Saunders  is  an  analyst  in  the  Data  Processing  Division  of  the  TRC  Service 
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analytical  groups.  His  duties  included  the  implementation  of  engine  tests,  test- 
data  analysis  and  performance  prediction,  and  programming  for  the  IBM  704. 

Educational  Background 
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APPENDIX  A 


SPECIFICATIONS  FOR  COMPUTATION  OF  FALLOUT  FROM 
A  TWO-DIMENSIONAL  CLOUD  IN  A  TWO-DIMENSIONAL  WIND  FIELD 


BLANK  PAGE 


APPENDIX  A.  SPECIFICATIONS  FOR  COMPUTATION  OF  FALLOUT  FROM 

A  TWO-DIMENSIONAL  CLOUD  IN  A  TWO-DIMENSIONAL  WIND 
FIELD 

A.  1.0  INTRODUCTION 

This  is  a  set  of  program  specifications  for  the  first  in  a  series  of  computer 
programs  to  be  used  to  evaluate  the  influence  of  local  wind  circulations  on  the 
deposition  of  close-in  fallout.  This  program  is  designed  to  treat  a  two-dimensional 
cloud  of  uniform-size  fallout  particles  initially  contained  within  the  first  few 
kilometers  above  the  earth’s  surface. 

The  program  was  written  in  FORTRAN  for  the  IBM  7094  digital  computer. 
Machine  time  requirements  may  be  estimated  from  the  following  result.  Three 
minutes  of  machine  time  were  required  to  deposit,  from  a  height  of  4  km,  eight 
wafers  (one  wafer  for  each  of  eight  particle  sizes).  Each  wafer  was  defined  by 
twenty-four  points  on  its  boundary. 

A. 2,0  INPUT 

Basic  input  will  be  from  cards  and  is  defined  as  follows; 


G 

Geostrophic  wind  speed,  cm  sec  1 

X.XXEXX 

A 

-X 

Ekman  spiral  parameter,  cm 

X.XXEXX 

XC 

Critical  separation  tolerances,  cm 

X.XXEXX 

YC 

Critical  separation  tolerances,  cm 

X.XXEXX 

XP(K) 

Coordinates  of  points  on  cloud  at  the  initial 
time,  cm 

X.XXXXEXX 

YP(K) 

Coordinates  of  points  on  cloud  at  the  initial 
time,  cm 

X.XXXXEXX 

KMX 

The  number  of  points  defining  the  cloud 

XMASS  - 
T0TMAS 

The  mass  of  fallout  material  within  the  cloud  at  the 
initial  time,  gms 

X.XXXXEXX 

DT 

The  time  step  to  be  used  in  computing  movement  of 
the  cloud,  initially  sec 

X.XXXXEXX 

T0UT(I) 

The  times  at  which  the  configuration  of  the  cloud 
is  to  be  outputted,  (I  <  50)  sec 

X.XXXXEXX 

VF 

The  fall  speed  of  the  cloud  particles  relative  to 
the  air  (negative  number),  cm  sec”* 

X.XXEXX 
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A. 3,0  COMPUTATIONS 


The  main  program  consists  of  a  series  of  subroutine  calling  statements, 
plus  a  scheme  for  condensing  output  from  the  subroutine  DPOSIT.  Once  deposition 
commences,  the  subroutine  DPOSIT  outputs  deposition  values  after  each  time  step. 

For  convenience  in  evaluating  these  results,  a  scheme  was  developed  and  incorporated 
into  the  main  program;  the  scheme  takes  the  outputs  from  all  the  time  steps  and 
orders  and  condenses  them  into  a  single  output  in  tabular  form  and  with  proper 
heading  identifications. 

AIL  computations  are  performed  in  the  subroutines  except  for  the  ordering 
and  condensing  of  the  output  from  subroutine  DPOSIT,  wnich  take  place  in  the  main 
program.  All  parameters  are  stored  in  common.  A  description  of  each  sub¬ 
routine  is  presented  below,  followed  by  the  calling  sequence  logic  (Fig.  A-l). 

A. 3.1  Subroutine  XINPW 

The  purpose  of  this  subroutine  is  to  compute  all  time -dependent  parameters 
such  as  U,  YT,  S,  and  XB  which  are  input  for  the  WIND  program.  These  param¬ 
eters  will  be  computed  at  increments  of  1800  seconds.  T  is  the  time  step  in 
seconds,  and  we  let  TD  be  the  time  of  day.  Whenever  T  =  1800  sec,  then  the  time 
of  day  (TD)  will  be  increased  by  0.5  and  the  four  parameters  mentioned  above  will 
be  recomputed. 

Output  of  this  routine  will  be  stored  in  common  for  use  as  input  to  the 
WIND  routine.  An  initial  TD  will  be  given  and  T  =  0  at  the  start.  Compute  the 
following  to  be  used  as  input: 

S  =  (-0.88(TD)2  +  26.4(TD)  -  166]  105 

U  -  [  -0.143(TD)2  +  4.086(TD)  -  22.571]  105 

If 

10  ^  TD  <  12,  go  to  A 

12  ^  TD  <  18,  go  to  B 

18  <  TD  ^  20,  go  to  C 

A.  XB  =  u 

YM  =  0.7  X  104 
YT  -  [0.4 (TD)  -  4]  10^ 
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B,  XB  =  ( -  1.33{TD)2  *  40(TD)  -  288]  105 
YM  -  0.7  X  104 

YT  -  l-0.044(TD)2  +  1.33(TD)  -  8.8]  IQ5 

C.  XB  =  0 

YM  =  0.7  X  104 

YT  =  [  -0.4(TD)  +  8]105 


A. 3. 2  Subroutine  WIND 

The  purpose  of  this  routine  is  to  compute,  using  analytical  formulas,  the 
u  and  w  components  of  the  wind  field  at  the  points  describing  the  slice.  The 
parameters  required  as  card  inputs  are  Z^5  Z^,  Z^,  A,  G,  SP(K),  and  YP(K). 
The  parameters  Z^,  U(t),  X(t),  and  (t)  are  computed  in  subroutine  XINPW  and, 
through  the  main  program,  are  available  for  use  in  the  subroutine.  At  this  point 
it  should  be  mentioned  that  the  notations  used  in  the  computer  program  and  in 
the  text  of  this  report  differ  in  some  instances.  These  differences  are  listed 
below. 


Text 

Program 

Text 

Progri 

zh 

YH 

cr 

SB 

z 

m 

YM 

U(t) 

U 

ZF 

YF 

u 

UBAR 

a 

A 

u 

U(K) 

G 

G 

w 

V(K) 

X 

XP(K) 

4> 

PHI 

X 

XB 

Z 

YP(K) 

a 

S 

At 

DT 

T 

YT 

Z 

YB 

The  output  of  this  subprogram  will  be  the  quantities  U(K)  and  V(K)  of  the 
wind  components  in  the  x  and  y  directions,  respectively,  at  point  K  where 
1  ;s  K  £  KMX. 
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Compute  the  following  to  be  used  as  input. 


SB  - 


YF 


YB  - 


YF 


-  YT 
4 

+  YT 


Let  K  =  1 . KMX 

If  YPK  25  YH,  go  to  A 
YPK  >  YH.  compute 

PHI  -  exp  [  - (XPK  -  XB)2/2S2]  and  store, 

HD  =  [(PHI)U]/[3(YT  ~  YM)2]  and  store, 

XD  =  (XPK  -  XB)/S2  and  store, 

VM  -  (XD)  (PHI)  U(YM  -  YH)/2  and  store, 

VT  =  VM  -  (XD)  (HD)  [ -  (YT  -  YM)3]  and  store, 
YM  -  YH  YT  -  YM 

XN  =  U(PHI)  { - - -  + - - - — )  and  store, 

-AYPK 

UBAR  =  G  [  1  -  (e  )  cos  (AYPK)]  and  go  to  B. 


A.  Compute 

UHAT  =  0 

-AYPK 

UBAR  -  G  [1  -  (e  )  cos  (AYPK)] 

VK  =  0 
UK  =  UBAR 

B.  If  YPK  25  YM,  go  to  1 

YPK  >  YM.  go  to  C 
1.  Compute 

UHAT  =  U ( PHI )  (YPK  -  YH)/(YM  -  YH) 

VK  =  (XD)  (PHI)  U(YPK  -  YH)2/2(YM  -  YH) 
UK  -  UBAR  +  UHAT 

C.  If  YPK  25  YT,  go  to  l 

\rPK  >  YT.  go  to  D 
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1.  Compute 

UHAT  -  U(PHI)  (YT  -  YPK)2/(YT  -  YM)2 

VK  =  VM  -  (XD)  (HD)  [(YT  -  YPK)3  -  (YT  -  YM)3] 
UK  -  UBAR  t  UHAT 


D.  If  YPK  s  YF,  go  to  1 

-AYPK 

YPK  ^  YT,  then  set  U(K)  =  L»  [  1  -  (e  )  cos  (AYPK)] 

and  evaluate  scheme  for  EY  where 


YPK 

EY  =  JYT  exp 


(Y  -  YB)2  ' 
2(SB)2  • 


dY 


1.  Compute 

XF  =  75J~(SB)  exp  <  "  <YPK  "  YB)2/[2(SB)2}} 

UHAT  =  (XN)  (XF) 

VK  =  VT  -  XN/lV2?(SB)] 


A.3.S  Subroutine  PMOVE 

This  subroutine  computes  the  displacement  of  the  points  describing  the  slice 
being  moved.  Input  is  the  coordinates  of  the  points,  the  wind  components  computed 
using  subroutine  WIND  at  the  points,  the  fall  velocity  of  the  particles  (VF),  and 
the  time  interval  DT. 

The  output  of  this  routine  is  a  new  set  of  coordinates  describing  the  slice 
after  time  step  DT,  which  are  input  to  the  following  subroutines:  ACCRCY, 

AREA,  and  DPOSIT. 


A. 3 .4  Subroutine  ACCRCY 

The  purpose  of  this  routine  is  to  check  the  separation  of  two  adjacent  points 
in  both  the  x  and  z  directions.  It  is  desired  that  two  consecutive  points  remain 
less  than  the  specified  distances  XC  and  YC  apart  in  the  x  and  z  directions 
respectively.  If  these  distances  are  exceeded,  the  subroutine  PMOVE  will  be 
repeated  using  DT  --  l/2(DT).  This  iteration  will  be  done  no  more  than  three 
times.  If  after  three  iterations  there  are  some  points  still  too  widely  spaced, 
then  the  last  results  will  be  printed-out  and  the  computations  for  that  particular 


fall  speed  will  be  terminated.  The  XC  and  YC  values  are  card  inputs.  Output 
from  this  routine  will  occur  only  if  the  separations  rerr  ain  too  large  after  the 
third  iteration.  The  separation  between  consecutive  points  K  and  K  +  1  is  given 
by: 

for  K  =  1,  (KMX  -  1) 

DX{K)  =  ABSF  [XPN  (K  +  1)  -  XPN(K)) 

DY(K)  =  ABSF  [XPN  (K  +  1)  -  YPN(K)] 
for  K  =  KMX 

DX(KMX)  =  ABSF  (XPN(KMX)  -  XPN<1)] 

DY(KMX)  =  ABSF  [XPN(KMX)  -  YPN(l)] 

XC  and  YC  will  be  input  in  the  present  program,  but  subsequently  will  be  computed 
for  input  to  this  subprogram . 

If  the  accuracy  criteria  above  are  satisfied,  set 
XP(K)  =  XPN  (K) 

YP(K)  =  YPN(K) 
for  K  =  1,  ....  KMX, 

and  set  T  =  T  +  DT  in  which  DT  is  the  value  used  in  the  PMOVE  program. 

A.3.5  Subroutine  AREA 

The  purpose  of  this  subroutine  is  to  compute  the  area  of  the  slice  and  the 
density  of  the  material  contained  within  the  slice  after  each  time  step.  Input  to 
this  routine  is  the  coordinates  of  the  points  after  each  time  step  and  the  mass 
of  material  contained  within  the  slice.  Output  is  density  which  is  required  as 
input  into  Subroutine  DPOSIT. 

Compute: 

K  =  1,  ....  (KMX  -  1) 

A(K)  =  1/2  *  [XP(K  +  1)  -  XP(K)j  *  [YP(K  +  1)  +  YP(K)J 


Then 


K  =  KMX 

A(KMX)  =  1/2  *  (XP(1)  -  XP(KMX)]  *  [YP(KMX)  +  YP(1)J. 
KMX 

"  Z  A(K)  . 

K=1 


AREA  = 


We  then  recompute  the  density 


DEN  -  XMASS/AREA. 

Note:  at  the  beginning  of  the  main  program  set  XMASS  =  T0TMAS 
A. 3. 6  Subprogram  DPOSIT 

This  subprogram  computes  the  amount,  if  any,  of  the  fallout  during  the  time 
interval,  DT,  assigns  it  to  the  x-coordinate  axis,  and  recomputes  the  area  of  the 
cloud  above  the  ground  together  with  the  modified  mass  within  the  cloud. 

The  first  step  in  the  program  is  to  check  the  Y  coordinates  of  the  KMX 
points.  If  none  of  these  is  less  than  zero,  one  may  exit  from  the  program  because 
no  deposition  has  occurred. 

Otherwise,  we  compute  as  follows: 

*  Let  K  =  1,  ....  (KMX  -  1) 

If  YP(K)  <  0,  go  to  A, 

YP(K)  2  0,  go  to  B. 

A.  If  YP(K  +  1)  >  0,  go  to  2 

YP(K  +  1)  =  0,  go  to  1 

YP(K  +  1)  <  0,  go  to  1 

1.  Area  =  -0.5  *  [XP(K  +  1)  -  XP(K>]  *  [YP(K  1)  +  YP(K)J 

Dep.  =  Den.  *  Area  [assign  to  interval:  XP(K)  to  XP(K  +  1)  and 

printout] 

XMASS  =  XMASS  -  Dep. 

XPN(K>  =  XP(K) 

YPN(K)  =  0 
G0  T0  * 


2.  Compute  X 


YP(K  +  1)  ♦  XP(K)  -  XP(K  +  1)  YP(K) 
YP(K  +  1)  -  YP(K) 


Area  =  -0.5  *  [X  -  XP(K)]  *  YP(K) 

Dep.  =  Den.  *  Area  [assign  to  interval  SP(K)  to  X  and  printout] 
XMASS  =  XMASS  -  Dep. 

XPN(K)  =  XP(K) 

XPN(K)  -  0 
G0  T0  * 
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B.  If  YP(K  +  1)  >  0,  go  to  1 

YP(K  +  1)  =  0,  go  to  2 

YP(K  +  1)  <  0,  go  to  3 


1  and  2  .  No  deposition 
Set  DEP(K)  *  0 

XPN(K)  =  XP(K) 
YPN(K)  =  YP(K) 
G0T0  * 


3. 


Compute  X  = 


YP(K  +  1)  *  XP(K)  -  XP(K  +  1)  YP(K) 
YP(K  +  1)  -  YP(K) 


Area  =  -0.5  *  [XP(K  +  1)  ~  X]  *  YP(K  +  1) 


Dep.  =  Den,  *  Area  [assign  to  interval  X  to  XP(K  +  1)  and  printout] 


XMASS  =  XMASS  -  Dep. 


Set  XPN(K)  =  XP(K) 
YPN(K)  =  YP(K) 
G0  T0  * 


For  K  =  KMX,  use  the  same  procedure  as  above  but  replace  (K  +  1)  index  by  1. 


A. 4.0  OUTPUT 

The  input  data  is  printed  out,  together  with  sufficient  comments  to  identify 
the  run.  The  final  output  of  the  run  will  be  a  tabulation  of  deposition,  with  its 
assigned  intervals  for  a  given  slice  and  fall  velocity. 


44 


All  Input 
read  In 


Fig.  A-l.  Calling  sequence  of  main  program. 
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APPENDIX  B 

FORTRAN  II  LISTING  FOR 
IBM  7090/94  FALLOUT  MAIN  PROGRAM 


THt.  FoLLU*lN<»  I  i>  A  FOKlfifcn  ii  Ll^I  iNb  F  on  loM  /Oxu/94  F<*ll*OuI 
MM  IN  FHObMfm  -  FAi_l.-OUT  19t>3 

U  I  MtNb  1  ON  A  I  <  l  boo  )  «(J  1  <  1  uoo  )  .  UL  l  <  i  buJ  )  «C  1  l  1  300  )  «o  1  <  1  300  >  * fc.  i  <  l  300  > 
0  I  MtNt  I  ON  XP<oOo>»YP<  30  O  )  .  XPN  (  30  0  >  •  YPN<  300  >«O(J00)»V!  30  O  ) 

Ul  MLNb  1  ON  UXI  300  >  «  LY  (  300  )  *  MR  (  300  )  •  Oc.P  <  300  >  t  X ( 300  >  »H<  3  > 

COMMON  XP« YP.XPN  •  YPN .u. V  <  OX  «  uY  •  MW  .  OC.N  .  KrtX  *  KM  i  »K  .  UP 

COMMON  XMAbb.  T  OTMAb.  T  *  o<  T  *  N«  VF  .  0. L V  .  oo  T 

COMMON  YH.YM.YT.VF»Xb«t»»A,Yfci..jet.X.H»XC.  YC  »ucP 

fit AO  333.0tA.XC. YC 

PRINT  OJt 

PRINT  333.OtA.XC.YC 
999  Rt AO  tt4t.LTt.Tu.MORc  »VF 
PRINT  tic  1 

d'ti\  FORMAT  (  3UH1  INPUT  VALutb  l  T  3  .  To  .  MORt  .  VF  ) 

PRINT  iitt.LTt.  tu.MORt.  VF: 
tti  FORMAT! 14. F?.t. 13. tit. t) 

Rt  AO  444.M.  (XP!K>  .*.=  1  >WI 
PRINT  443 

PRINT  444 .M* (XP l X) «K»l .M ) 

RtAO  444.  M.(YP(N).M1.M) 

PRINT  44i> 

PRINT  444. M. ( YP t X ) «K*1 .M) 

RtAO  33b.  TQTMAb.OT. YH.YF ,XMX 
PRINT  334 

PRINT  3b3.TOTMAt.OT . YH.YF .KMX 
0O9UUKA® I . lboo 
A  1  (KA) so. 

31 (KA) so. 

Otl (KA»so. 

900  CUNT  1 NUt 
K0NUM=0 
MA»KMX 
MTsMA 
KMl«KMX-l 
XMAtbsTOTMAS 
NT  =0 
T  *0  • 

t  nt*nt+i 

4  IF(NT-LTb>3«3* lt»U 
3  CALL  XlNP«(TO> 

6  CALL  Wl NO 

7  N*0 

OUT st. *uT 

8  DOT «  »3*oOT 
CALL  PMOVE 
CALL  ACCRCY 
1F(N-1 >70.3.9 

9  IF(N-3)b.t«»6.t66 
dt>b  PRINT  500.  T. NT.  TO 

DOt  e»«K  ■  1  •  KMX  .  3 

LOB  PRINT  301 .XPNtKJ .YPN(X) .XPNlN+1 > .YPNIX+i J .XPNIKtt) .YPNlUtc > 
lF(MORt-l >130.999.999 
70  KT«T+OUT 
CALL  ARttA 
D0100K-1 *KMX 
IF ( YP( K> ) 105.100. 10O 
100  CONTINOL 
GO  TO  IOC 
105  DO103K-I.KMX 
OtP(K) so. 
loo  X(K>»U. 


CALL  UPUblT 
OGIOOX* 1 • KMX 
XP(K>*XPN(K  ) 

1 09  YP(X > » VPNCX i 

XPIKMX+i >*XP( 1 ) 

YPCKMX+i  >  « VP  l  l  J 
NK*0 

XUNUM*KUNUM+1 
00646*= l • KMX 
1F(DEP<X)  ) 664 *646 *634 
«>S4  I F  <  XDNUM-  I  >  656  .6 66*657 
6b?  IF  (A  1 1  X  )  -XP(  X  >  )  650*  633*65-» 

63b  1F(6I(KJ-XP<X+1  J  163U*63Y*63-> 

633  IF(A1 (K> 1636*666*633 
638  NK*NX+1 
MT*MAf NX 
Ai (MT ) *XP(K) 
ailMT) *XP(<*1 ) 

ot  1  (  MT  )  *UE  1  (MT  )  +DtP  (  K  1  ✓ AtJSF  <  XP  ( X  1  -XP  (  X+  1  )  ) 
&0  TO  646 
666  All  K)*XP(K> 

QIC  K)*XP(K>1) 

669  Dfc 1  (  X  ) *U£ I <  X 1 +OtP<N 1/A3SF ( XP (X > -XP( X+ i > ) 
646  CONTINUE. 

MA*MT 

IF  <  MT  — 1 3uO 1 130* 160* I 6U 
1 30  D0140X*1*KMX 

IF  I YP  (X) 1 13u* 140* 130 

150  T«T+OuT 
IF(T-iauo. 12* 151 « 151 

151  TO*TO+.6 
T*U* 

1 F  C  TO— HU*  >  £  *2  * 1 60 
140  CONTINUE 

D0675K*! *MT 

1F(A1(X>-b1(N) 1675* 676 »6?7 
677  A 1 A* A 1 ( X ) 

B 16*6 1 C  X ) 

Al (K)*616 
B1 (X ) *A 1 A 
675  CONTINUE 
M*2 

7U3  ZMAX=A1(M-1) 

ZMAX1 *81 (M-l  ) 

ZMAX3*0fc 1 ( H-  1  ) 

00701  I *M.MT 

IF(ZMAX-Al < I J ) 701 »7wl *704 
702  ZMAX*A 1(1) 

ZMAX 1 *6 1 4  I  ' 

ZMAXO*t)fcl  (  1  ) 

K 1  *  I 

701  CONTINUfc 

SAVA  1 *A 1 (M-l  > 

SAV6 1 “6 1 <  M— 1  1 
SAVOfc l*ufcl (M-l ) 

AI (M- l > *Z MAX 
B 1 (M-l 1 *ZMAXl 
DEI <M-l >*ZMAX3 
Ai  (XI  )*3AVA1 
Q  1  (XI ) *5A VO 1 
otl (XI ) *SAVUfcl 
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M  =  M+  1 

IF  (M-MT  )7oJ«7oJ,7o4 

714  C 1  (  1  >  =  A 1  (  1  ) 

Q 1 ( 1 ) =6 1 ( 1 ) 

El  <  1  > =  OE 1(1) 

1  I  =1 

JJ  =  2 
MM  =2 
CC  =  2 

715  I F ( A  1 <  J  J ) -U 1(1I))71I*71I,71<£ 

712  Cl  (EL.)  =01(11) 

01  (LC)*Al  <  JJ) 

El  <LL)«0, 

Cl ( LLt 1 IxAl (JJl 
01 ( CC+ 1 1*01 <  JJ ) 

El  (CC+l  )«OEl  (JO) 

CL*CC+  1 
GO  TO  706 
71  1  Cl  (  LC  )  *  A  1  (  JJ) 

0 1 ( LC ) *b 1 ( J J ) 

El (LL) =uE 1 (JO) 

706  11=11+1 
JJ*JJ+1 
LL=Cl_+  1 

IF  (MM-MT  >713*714*714 

713  MM  =  MM  + 1 
GO  TO  716 

714  l_MX  =  EC  —  1 
CMX2=LC-2 
DOSOol*l  ,CMX2 

IF (Cl ( 1  +  1 ) -Cl (  1 ) >  80 1 .602*803 
SUE  A l ( 1 >«C1 ( 1 > 

Bill) =C 1(1+1) 

OEM  1  )  -0. 

601  GO  TO  600 
Qo3  Del < i ) =o. 

A 1  ( 1  ) =C  1  (  I  ) 

61 < I )=Ci < 1+1 ) 

0U6U4KK= 1 , 1 

IF <01  IKK) -Cl ( 1  +  1  ) ) 606 *605.603 
605  DEI ( I >*UE1 < 1 )+El <KK)*(C1 ( 1+1 )— Cl ( l ) ) 

GO  TO  804 

806  IF ( O 1 ( KX ) — C 1(1)) 607  «  807  «  606 
6o7  0E1 < 1 ) =oE 1 ( 1 ) 

GO  TO  604 

808  OE  1  <  I  )  =uE  1  <  1  )  +c  l  <  K.K  )  *  (  u  1  (  KK.  )  -o  1  (  1  >  ) 

604  CONTINUE 
600  CONTINUE 

A 1 ( CMX ) *C 1 ( CMX ) 

61 1 CMX ) «u 1 (CMX ) 

OE 1  ( CMX  >  *E 1  (LMX) *(Ul  (CMX ) -C 1  (CMX  )  ) 

0061 OME* 1 .CMXc 

IF (01 (ME) -Cl (CMX) ) 61 0*61 0,611 
811  OE 1 <  CMX ) =OE 1 ( CMX )  +  E l  <  ME ) * ( O 1  ( Mt ) -C 1 <  CMX  )  ) 

610  CONTINUE 

PRINT  6oj ,KT  * TD«NT 
PRINT  6cO*VF 
D0925K*! • CMX 

925  PRINT  6JO.K.OE1  (K)  *A1  (K)  *61  (n) 

IF ( MORE- 1 ) 160.999,999 

332  FORMAT ( 46H  THE  FOCCOalNO  OATa  are  INPUT  VAlolj  g.A.aC.YC) 
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333  FORMAT <4 El  2. 2) 

444  FORMAT  !  13/(6c.12«E)  ) 

443  FORMAT  <4^H  T  Mfc  FOELQwiNG  DATA  ARL-XP-InPgT  VALUES  > 

445  FORMAT (4UH  THE  FOLLOWING  OATA  ARE-YP-INPuT  VALUES) 

500  FORMAT  <  S5H1  THIRL  ITERATION  HtACHtu.  ACCURACY  GREATER  THAN  XC.YC 

123H  TIME  IN  SECONUS  s  FU  •  0/  A  X  .  JHXPN  i  YX  ,  JH  YPN  t  YX  *  JHXPN  i  VX  .  Jri> 

2N.9X  .3HXPN.  OX*  UHYPN/  i  7H  NuMctcW  uF  PASO  **  14/lJH  HOUR  OF  uAYaiF7*£ 

501  FORMAT ( 6E 14*51 

554  FORMAT! D9H  THE  FOLLOWING  DATA  ARE  INPUT  VAloEe  TuTMAl * u T * YH » YF * 

1  KMX  ) 

555  FORMAT ! 4E 12.2. 13 3 

600  FORMAT! IhI . JX. IoHTIME  IN  SeCONuS =. i 6/ 1 3H  HOUR  OF  DAY * »F 7 *2/ 1 UM  Pa 
IS  NUMBER* « 143 

601  FORMAT! I3.7Ei2.s3 

620  FORMAT i X/46H  SUMMATION  OF  OEPOSIT  WITH  ITS  PROPER  lNTERVALXOH  V 

1 .E12.2/4H  K.6X.7HUEPOSIT. 14X.2HA1 , 16X»2Htil ) 

630  FORMAT! 14. SEItt.S) 

160  CALL  EXIT 
END 
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OULjNOOl  INL  X  jhtf'nUul 

U  1  MtNj  i  ON  XP  <  Ju  j  *  .  *  P ( JOu i  *  XPIv  (  JUUl  .  VW-I  JuO  I  «  ut  JOtll  <VlJiiU) 
U  I  MtNj  i  UN  OX i OUu )  lurUvu)  .AM(  Jlyu  J  iWi.K'i  JUUI  iXU^viirti  J) 
L.UMHUN  XP  «  YP  •  XPN  *  yPUiUi  V  tuXtuV  lAHiwLti.XMx  tUMl  i*«uP 
COMMON  XMAjj  ilU!  Mi  -  iT  •  u  I  «  N  «  Vt-  •  C  •  O  Y  •  uu  T 
C  OMMON  YH  *  YM  i  Y  I  •  Yf  *Yd*-jU«XiHiXC*  YL 

It-  I  TO-  lu>  >  l‘jb<  1  JXi  J  j.; 

13£  TUi»*TU**x 

i»»  <  t  —  •  bo*  f  ub  i  •*  j  *  i  •  Ot.  o 

OP*  <  i  -•  I  43*TuO  >  +  4»  UU6*  I  o -£.*-•  o  /  I  }  *  1  ■ix.< 

IF  (  TO-  ix»  J  1  jA  ,  i  JA 

13*  IF  |  TU-  lb.  1  106.  lot,»  1  JO 
13b  IFtTo-xu. Mb/. lu7, l^b 
IbS  X0*0*0 

T«»  I  .4  J  *  1  .  OL.O 

YT  «  <  •*•  i  0-«*  >  *  )  •  OL.O 

OU  T  U  1  Ob 

1 06  xb*  <  —  l • vj*tu-+Au>o* I u*xob • t * i • otb 

YM*  <  .4 ) *  1  .ULO 

YT*(-«U*4*fL)j+  1  .  JJ*!U-U»U1  *1  »OLS 
OO  T  O  1  Ob 

107  xe=  o»o 

YM  =  <  *4 ) * l  . UCO 
YT=<— •4*TO+U«  )  *  1  .OLO 
ibb  PbTUPN 
tNU 
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3UO ROUT  i  Ni_  *  i  NL 

U  i  MLN;»  1  UN  J'PUtvlirFUUv!.  XHNl  Jv(/>  •  YPM  JOO  )  «‘vl  Jgu  )  *  V  (  _JU  U  ) 

DI  Mfc!  ii>  i  UN  UXUiJU  )tUV<  300  l.Afil  ->Qu  >  tUt-P  <  J'JU  )  tXt  JQQ  i  itllJ) 

COMMON  XP, YP . XPN « YPN. U.V. OX. UY  , A« » ULN » KMX  *  KM  1  , K i UP 

COMMON  XMAh>i>  *  T  OTMAu  1T1  DT  « N t VF • G * LY ♦ UOT 

COMMON  YN  »  YM«  Yt  .  YF  •  XB*  3  •  A  ,  YO  .  ati  ,X  «H.XU  ,  YC  t  UtP 

P1*J.  1415927 

3d  «  (  YF  -  Y  T  )  /  4  , 

YB*<  YF  +  YT  >/2» 

DO  40  Ns l , KMX 
AYPK=A*YP<K ) 

ei  *i ./expfi aypk j 

B2=CG3F  <  A yPK  ) 

IF  1  YP(X  F-yh  )  20  .20*  1  J 

13  pi = c xp< n j -xa i **a 

P2*2.*S**2 

PHi*l «/LXPF <Pi/P2> 

HJ  *PH|*UP 
YTM* YT-YM 
H2=3.*YlM**2 
HO*Hl/H2 

XU*IXP(K)-XttJ/b**<i 

YMM=YM-YH 

VM*XD*Hi*YMH/2. 

VT=VM-X0*<-YTM**J  )  *HO 
XNIsYMH/i.+yTM/j. 

XN*Hi*XNI 
UBAR«G*( 1 .-01*02) 

GO  TO  23 
20  UHATso. 

U8AR*G* ( 1 .-01*02 > 

V<K) *G. 

U<K)*UOAfi 
GO  TO  40 

23  IFtYP<K)-YM)2t».25.30 
25  UHAT*H1*<yP<K)-YH)/YMH 

V<K)  »XU*MI*(  <  YP(N)-YM>**2  )/<«;  •  *  ymh  ) 

UIK  )  =UBAR  +  UMAT 
GO  TO  40 

30  IF  <  YP (X)-yT ) 3j • 33  «  35 

33  UHAT*H1*I ( YT-YP(K) ) **2  )/<YTM**2  ) 

YTKM-I YT-YP<K) >**3  -YTM**3 

V  <K I «VM-XD*H0*YTKM 
U<K ) *UOAR+UHAT 
GO  TO  40 

35  IF«YP<K>-yF)38»30»44 
44  V(K>*0. 

U<K>  *G*  t 1 .-01*02 ) 

GO  TO  40 

38  YKOS*< YP<K)-Yo)**2 
T 503 ■  2  *  *  <  80**2  I 
XF i * 3QR T F  <  2 • *P I >*SO 
XF2*  J ./XFi 

XF*xF2*  (  1  ./EXPF<  VKos/T bUb  )  ) 

UHAT»XN*XF 
caul  sintey 

VIK ) *VT-<  XN/XF l l*EY*XO 
U<K  ) aUOAR+UHAT 
40  CONTINUE 
RE TOWN 
ENU 


54 


J  i  nc.  In  ^ 

U  t  Mt .Ni  2  uN  (  JUui  t¥^i  JUj  )  ,  xHN  t  JwUi  tYPM  OuU  )  lUUUDi  •  V ( 30  O  > 

D  I  MENS  :  UN  UX ( jQU )  *l/V(  JUU)  .  AK  (  ^ow  )  tut-p  t  300  i  *Xt  OQO  i  t  H  t  J  ) 
COMMON  XP.YPi  XPN  *YPNiU«V<UX»uY  •  **R  .ut.NiK.MX  «  KM  1  .  N  •  UP 
COMMON  XMAss  ,  f  OT  MA^  *  T  «  U  I  *Ni  Yr  il)  .tY.UwJl 
COMMON  YH.YM.yT.YFiXb.J.A.Yc.jO.X.MiXt, YC.utP 
3  INK  =  -J<_« 

YLL - «  T 

1  NK  ai  1  NK/-.E  . 

2  1 NK  -  t  YP(K)-YT  )  /  3  I  NK 
SUM=0. 

R  =  U. 

CONsi;.*  <  3u*3o) 
to  l  a  I  =  1 »  IKK 
OO  1  t>o  =  i «  3 

TLMP=( t <Y£t+R)~YB)**3  3/CON 
H<  U  >  =  1 ./£XPF «  TtMP > 

16  R=«+XINK 
•m!  S  U  • 

SUM-3UM+f7.tNK/J.  )  *  t  rt  t  1  )  +  4  a  *rt  (  A  )  +  H  <  3  )  J 
15  YtLsYtt  +  c.tilINK 
E  Y*SOM 
RETURN 
END 


SUBROUTINE  PMOVE 

O  I  MENS  ION  XPiOOO)  tYP<300)  .  XPNl  300)  .YPN1  JOQ  >  .u<3UO>  *VUO  U) 

DIMENSION  DXOOU) • DY  < 300 ) «AR(30U> »UEP(SOOJ  .X(JOO)  .H<  J) 

COMMON  XP . YP . XPN . YPN  «  U  «  V  «  DX «  UY  «  AR • UtN  « KMX .KMl.K.oP 

COMMON  XMAsS  «  T  OTMASt  T  «  UT  »N«  VF  *  G«  EY iUoT 

COMMON  YH • YM.YT  *  YF  »Xto»3« A , YO . SB » X t M »  XC , YC « D£P 

DU  4^sl,KNX 

FKT  =- ( V  t K ) +VF ) *UOT 

IF  t  YP(K)-FKT  138*06.  Jt> 

36  XPN(K) =XP(K)+u<N> *UUT 

YPNIM  =YP(K)  +  {  V<K.J+VF  J  *UOT 
GO  TO  4£ 

3B  UTP=-YP<K>/(V(K)+VF) 

TAU*DOT-DTP 

XPNCK) »XP(K)+UtK»*OTP 

YPN(K) =VF#TAU 

4£  CONTINUE 
RETURN 
END 
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UUOwUUT  i  Nfc  A4~i_RL.  ¥ 

D  i  Mi_N-a  i  s-N  XP(  juu)  •  YP  (  Juu  )  «  XPN  4  juul  *  VPN  (  JUU  I  tVl  Jutl)  »  V ( 30  0 ) 
L)  i  MENS  1  UN  UX(JOO)  Il/Vt  JUuNAWt  JOW(t¥tP'l  J¥¥l  •  X  !  J06  *  *  J  ) 
COMMON  XP  .  Y1J  iXPN  *  VPN  i^iVi  D*  <  uY  .  m  R  ,  uuN  .KMX  *  Km  i  <  K.  t  UP 
COMMON  XttASS «TOTMAu.T*Ut*N*VP  *  C  .  C.Y  •  UuT 
COMMON  VH»YM<YT.VK*Xe*a*A(Yo«JB.X.M.X4..YC*  Dt.P 
DO  4t*K  = 1  .KM1 

OX  (K  )  *At>SF  4  XPNIX  +  1  )  -XPNiK  1  ) 

1 F  lXC~DX(Xl  )  t><_  ♦  6<^  «Ao 
46  DY (K > =Ao5F ( YPN(K+ 1 ) ~YPN(K ) ) 

JF(YC-UY(K>  IbSiOt.'U) 

4  B  CONTI NUt 

OX  (KMX  }  =  Abi,F  l  XPN (KMX ) -XPN (  1  )  ) 

1 F  (  XC-DX  (  KMX  )  )  6E  .  3*-  *  50 
50  DY(KMX) =AbbF( YPN(KMX J-YPNI 1 > ) 

IF( YC-DYIKMX)  )63i6<i.54 
62  N=N+1 
RETURN 
04  CONTINUE 
N-0 

DO  55K=1.KMX 
XP(K ) = XPNIK ) 

55  YP(K»aYPN(K) 

RtTURN 

END 


SUBROUTINE  AREA 

DIMENSION  XP(3GQ) .YP(300> ,XPN( 300) .YPN (300) ,o( 300) .V(3QG> 

DIMENSION  ox (300) »DY(300) , AR ( 300 ) . OEP ( 300 ) »X(300> .H(3) 

COMMON  XP.YP. XPN . YPN  «  U . V  «  DX . UY . AR . DEN . KMX . KM 1 . K . UP 

COMMON  XMASS  »  T OTMAsi T  *  UT . N . VF . G . E Y » DOT 

COMMON  YH.  YM»  Y  T  .  YF  .XB.3.  A,  Yt>  .  ob.X.H.  XC  .  YC.UEP 

OO  75K* 1 . KM 1 

DXP*XP(K+l )-XP(K) 

OYP*YP(K+i )+YP(K) 

75  ARik) =.5*0XP#DYP 
OMXP*XP (  l > -XP  <  KMX ) 

DMYP=YP(KMX)+YP( I ) 

AR (KMX ) =.5*0MXP*0MYP 
AREA«U. 

DO  76K=l.KMX 

76  AREA*AREA+AR(K ) 

OEN=XMA3S/ ( -AREA ) 

RETURN 

ENO 
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GUOPOwTlNt  uPO-sIT 

Q  J  MENU  i  WK  APi  JUOi  i  YH{  Juu  I  JwiJ  J  *  YP^i  JOU  )  l  JVW  i  *  V  i  JUU  ) 

U  1  MLNS  1  UN  UX  (  JUU)  »UY(  JOU  )  *  *»N  <  ->Uu  )  *  utH  (  -SUM  >  «*(  JO'S  )  ♦  H  (  u  S 
COMMON  XP  *  YP(  XPN  *  YPNtutV  «  OX  *  U  T  .  MH  »  u{.  N  .  K  i-l*  *  Sh  i  *  X  *  UP 
COMMON  XWAijtToT  MA  u  «  T  ♦  <_>  T  «N«VF«o»t.Y*uiwT 
COMMON  YM,YM.YT.YF.XtS»0»A.YC3tUD.X«M.XU.YC«  UtP 
XPIKMX+ 1 ) =  XP< 1 > 

YP  <  KHX+ 1  ) sYPI  i  } 

DO  i  25K  2  1  «  KMX 
KK*K+ 1 

joe  i r «  yp(k)  j  i  io«  i  its* i  is 

110  IF<  YPIKM  )  1  l  1  «  1  1  1  *  lt£ 

1  l  1  AHfc.1  »-.S*(XP(AK  1  -XP<K)  )*(yP(NM  +  YP(M) 

O tPLK ) *utN*ARt 1 
XMAiSaXMA5i-UtP  i  K ) 

XPN(K) sXP(K ) 

YPN1K.J  *U. 

GO  TO  125 

1  1  2  YXXY=YP INK > *XP<K )-XP (KtO  *YP  <x ) 

XCK 3  =YXXY/<  YP IKK ) -YP (K  )  > 

AWE1  =-.5*<X{K)-XP<K)  )*YP<K  ) 

0EP<K)sUEN*AREl 
XMASE=  XMA5E-OEP  (  K ) 

XPN(K) *XP<K J 
YPNtKl s0. 

GO  TO  1a5 

lie  if  ( yp(kk) > i2o. i ie. i ib 
i ie  dep(K)so. 

XPN ( K ) *  XP ( K ) 

YPN(K)»YP<K) 

GO  TO  126 

1 2 0  YXXY* YP ( KK  > *XP l K ) -XP ( KK 3 *YP ( K ) 

X(K)*YXXY/t YP(KK J-YPtK 1 ) 

ARE 1  *—  • 6*  <  XP ( KK | —X ( K ) )*YPCKK) 

DEP(K)«UEN*AR£1 
XMAS6*  XMA66-0EP <  K ) 

XPN ( K ) s  XP IK) 

YPN (K)*YP(K) 

125  CONTINUE 
RETURN 
ENU 
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appendix  c 

THREE-DIMENSIONAL  FALLOUT 
COMPUTATION  SCHEME 


I  sm  K  i  L>f*U 


:> pr \  f)|  \  {. '  t  hr K  t  LH Si  t 

The  stem  and  cloud  is  origin 


!!%  \  !  ?  A  1  l.nr  i  J  ATIUNA1.  SHiKMi. 

all\  sulwii  v  idl'd  into  discs  and  further  sub¬ 


divided  into  wafers.  As  was  stated  m  the  text  of  this  report,  when  the  center 
of  the  wafer  reaches  an  altitude  below  which  local  circulations  arc  permitted, 
we  further  subdivide  the  w'afer  into  smaller  volume  elements.  Presented  below 
is  the  method  by  which  this  further  subdivision  is  achieved. 

Given  a  cylinder  C  -  BXI  where  B  is  a  disc  of  radius  R,  and  I  -  (0,  d), 


z 


C  is  to  be  partitioned  into  equal  volume  elements  each  of  which  has  a  “center  of 
concentration”  whose  coordinates  are  to  be  determined.  The  partitioning  is  to  be 
done  by  taking  m  +  1  cross-sections  at  right  angles  to  the  z-axis  to  give  m  layers 
of  height  d/m,  and  then  by  distributing  the  centers  of  the  volume  elements  in  each 
layer  as  follows. 


y 


If  we  take  concentric  circles  of  radius  (R/n)j,  j  =  1,  n,  we  can  decompose  any 

cross-section,  A,  into  the  sum  of  n  annuli,  . A^  (the  center  degenerate  annulus 

A^  ueing  the  (R/n)-disc).  The  area  of  the  jth  annulus  is  given  by 
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R 


Area  (A  )  t  (  ~  j ) 
J  » 


R 


(J 


,R,  .,2 

=  ^nj  (J 


R 


2 


-  *(“)  (2j  -  1) 


n  R  2  n  R  2 

Note  that  £  Area  (A,>  -  jr(— )  ]T  (2j  -  1)  -  ir(  —  )  (1  +  3  +  ...  +  2n-lj 

j=l  J  n  j= 1  B 

9 

R  n  2 

=  t( — )  T(2n)  =  itR  ,  the  total  cross -sectional  area.  Area  (A), 

n 


Then 


R 


Area  (A.)  jt(-) 

_ L  _  -.‘.D-..-  i2j...7..D  _  2j _ r _ i 


Area  (A) 


<  1 


ttR  n 

gives  the  fraction  of  the  total  area  taken  up  by  A 


J 


If  in  each  layer  there  are  required  to  be  N  centers  of  elemental  volume 


ivouli  be  desirable  to  have  [{2j  -  l)/n  }  •  N  canter  points  in  the  jth  annulus 


A^ ,  j  1 ,  • » . ,  n  , 

Then 


N.  = 
J 


ZLjzJ, 

2 


N 


n 


2  2 

must  be  an  integer  for  each  j  .  Then  N/n  must  be  an  integer,  or  N  -  Kn 


for  K  an  integer.  So  given  the  overall  number  of  points  N  ,  we  must  find  suit; 
K  and  n  such  that  the  preceding  relation  holds. 


In  each  A.,  j  /  1 ,  it  is  desirable  to  distribute  the  “center”  points  or 


mid-circle  of  the  annulus  at  equal  distances  apart.  These  mid-circles  have  i 


(R/n)  (j  -  1/2),  j  =  2,  n.  In  the  jth  annulus,  the  ISb  points  are  27r/?b  radi 


apart.  If  one  point  is  located  on  the  x-axis  (so  that  0  =  0)  in  polar  coordinati 


the  “center  points”  in  A^  are  at 


(P.©)-  = 


*  «  - 1/2)’ 

J 


P  =  0,  1 . N 


j  =  1,  2,  ...,  n  . 
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Fur  i 
S 

1 


1  (the  (.enter  discs,  v>v  put  one  point  at  the  origin  UK  U»  and  the  remaining 
1  {Knot 8  at 


R  2  sr 


p  -  0, 


so  that  they  lit*  on  the  circle  of  radius  1/2  (R/n)  . 


